
function [traces,TT]=PassiveFW(model,source,S_R_Geo,flag);

% written by Zhang Kai especially for passive surface wave modelling
% IGGE 2021/9/30
% Parallel calculations over sources

%**************************************************************************
% The main idea on modeling seismic full wavefield is based on
% transversely isotropic layered medium.
% So the source-receiver geometry (distance and azimuth) is the only variable that
% contributes to wavefield.

% We calculate each wavefield generated by individual sources
% and then convert cylindrical coordinates to Cartesian coordinates. Finally,we stack them together
% according to their random sources' excitation time.
% Caution should be taken: Each shot gather must have enough time to completely record the seismic
% signal within the corresponding source-receiver distance.

% This version is suitalbe to both the elastic and viscoelastic media.
%**************************************************************************
% input parameter
% delka_t = duration of the whole synthetic waveform  [s]
% source_number = ambient sources number
% dt = sampling frequency - in time [s]
% fmax=max. frequency to be calculated
% S_R_Geo= Source-Receiver geometry: Sources distribution and Receivers distribution

% dt and fmax is the key factors that affect computation efficienty.
% dt=0.002 is accurate enough for simulating near wavefields, while
% dt=0.01 is a better choice for >100 m wavefields.
%**************************************************************************
% output parameter
% traces.z: ambient noise,vertical component
% traces.x: ambient noise,x-axis component
% traces.y: ambient noise,y-axis component
% TT.offset: distance from all sources to all receivers
% TT.r: seismic record at all receivers generated by each source,radial component
% TT.t: seismic record at all receivers generated by each source,tangential component
% TT.z: seismic record at all receivers generated by each source,vertical component
% TT.t_one: time series corresponding to TT_z
% traces.t_total: time series corresponding to traces_z
% TT.fm: main frequencies of all sources
%**************************************************************************
% Material properties (size of arrays = number of layers)
%**************************************************************************
Cs=model.vs;     % Shear wave velocity [m/s]
Cp=model.vp;
Rho=model.dns;    % Mass density of layers (kg/m^3)
H=model.thk;      % depth of layers [m]
Damp_p=model.Damp_p;%1./(2*Qp); % damping
Damp_s=model.Damp_s;%1./(2*Qs);
%**************************************************************************
%% source and receiver parameters
source_number=source.number; % ambient sources number
delka_t=source.recordlength; % record time(s)
dt=source.dt; %  time step
fmax=source.maxfre;  %% max frequency to be calculated
fm1=source.fmin;  %% min. Dominant Frequency of Ricker wavelet
fm2=source.fmax;  %% max. Dominant Frequency of Ricker wavelet
%%%%%%%%%%%%  set point loads   %%%%%%%%%%%%%%%%%%%%%
Fx=source.Fx; % point load in x-axis
Fy=source.Fy; % point load in y-axis
Fz=source.Fz; % point load in z-axis,z-axis is positive in the upward direction

fm=rand(1,source_number)*(fm2-fm1)+fm1*ones(1,source_number);  %% Random distribution of Wavelet Dominant Frequency
delay_t=4/fm1;             %% excitation time of delta function
%**************************************************************************
%%%%%%%%%%%%  set ambient source-receiver geometry    %%%%%%%%%%%%%%%%%%%%%
sourceX=S_R_Geo.sx;
sourceY=S_R_Geo.sy;
receiverX=S_R_Geo.rx;
receiverY=S_R_Geo.ry;
ambient_geometry_plot(S_R_Geo,1); %% plot

N_re=length(receiverX);
y_pos_all=sqrt((sourceX-receiverX').^2+(sourceY-receiverY').^2);   %% compute offset
%**************************************************************************
y_max=ceil(max(max(y_pos_all))); %% max. range of Y coordinate
tmax=2*y_max/min(Cs)/0.862+2/fm1+delay_t;  %% max. duration of one record
srt22=ceil((delka_t-tmax)*rand(source_number,1)/dt);  %% sources' excitation time
t=0:dt:tmax;
N=length(t);
N2=0;%fix(y_max/min(Cs)/0.862/dt); %% mute the tail of the record
if mod(N,2)>0   %% make N even
    N=N+1;
end

N_shift=fix(delay_t/dt);
t_total=0:dt:delka_t;
traces_x=zeros(length(t_total),N_re);
traces_y=zeros(length(t_total),N_re);
traces_z=zeros(length(t_total),N_re);
TT_r=zeros(N-N_shift-N2,N_re,source_number);
TT_t=zeros(N-N_shift-N2,N_re,source_number);
TT_z=zeros(N-N_shift-N2,N_re,source_number);
t_one=0:dt:dt*(N-N_shift-N2-1);
%% calculation of green function in f-k domain, only once
UK=Green_HV_delta_one_par(model,N,dt,y_max,delay_t,fmax,flag);

offset=zeros(source_number,N_re);

%% start of all receivers due to each point source and then stack
numIterations = source_number;
cpu_N=feature('numCores');
str_kai=strcat(['Parallel computation using ', num2str(cpu_N)], ' cpu cores, Please wait...');
hbar = parfor_progressbar(numIterations,str_kai); %create the progress bar

parfor ii=1:source_number
    ut1_v=0;ut2_v=0;ut3_v=0;
    ut1_h=0;ut2_h=0;ut3_h=0;
    
    y_pos= y_pos_all(:,ii);
    y_pos=y_pos';
    offset(ii,:)=y_pos;
    [Azi,~]=cart2pol(receiverX-sourceX(ii), receiverY-sourceY(ii));   %% compute azimuth
    
    %% do the Hankel transform from f-k domain to f-space domain using discrete wavenumber method
    [Wavef,tt]=Green_HV_delta_two(UK,N,dt,y_pos,fmax,flag);
    %*****************************************************************************
    if any(strcmp(flag,{'v','V','hv','HV','Hv','hV'}))
        % z-axis point load
        % Vertical component
        ut_temp=conv_wavelet(tt,Wavef(:,:,5),fm(1,ii)); %% calculation of all receivers due to each point source
        ut3_v=ut_temp(N_shift+1:end-N2,:);
        ut3_v=mute_kai(ut3_v,1.5*max(Cp),y_pos,dt); %% mute events faster than 1.5*max(Cp)
        
        % Radial compenent
        ut_temp=conv_wavelet(tt,Wavef(:,:,4),fm(1,ii)); %% calculation of all receivers due to each point source
        ut1_v=ut_temp(N_shift+1:end-N2,:);
        ut1_v=mute_kai(ut1_v,1.5*max(Cp),y_pos,dt); %% mute events faster than 1.5*max(Cp)
        
        % Tangential compenent
        ut2_v=zeros(size(ut1_v));
    end
    %*****************************************************************************
    % xy-axis point load
    if any(strcmp(flag,{'H','h','hv','HV','Hv','hV'}))
        % Vertical component
        ut_temp=conv_wavelet(tt,Wavef(:,:,3),fm(1,ii)); %% calculation of all receivers due to each point source
        ut3_h=ut_temp(N_shift+1:end-N2,:);
        ut3_h=mute_kai(ut3_h,1.5*max(Cp),y_pos,dt); %% mute events faster than 1.5*max(Cp)
        
        % Radial compenent
        ut_temp=conv_wavelet(tt,Wavef(:,:,1),fm(1,ii)); %% calculation of all receivers due to each point source
        ut1_h=ut_temp(N_shift+1:end-N2,:);
        ut1_h=mute_kai(ut1_h,1.5*max(Cp),y_pos,dt); %% mute events faster than 1.5*max(Cp)
        
        % Tangential compenent
        ut_temp=conv_wavelet(tt,Wavef(:,:,2),fm(1,ii)); %% calculation of all receivers due to each point source
        ut2_h=ut_temp(N_shift+1:end-N2,:);
        ut2_h=mute_kai(ut2_h,1.5*max(Cp),y_pos,dt); %% mute events faster than 1.5*max(Cp)
    end
    % sum displacements over ur,ut,uz in cylindrical coordinates
    
    ut1=(Fx(ii)*cos(Azi)+Fy(ii)*sin(Azi)).*ut1_h+Fz(ii)*ut1_v;
    ut2=(Fy(ii)*cos(Azi)-Fx(ii)*sin(Azi)).*ut2_h+Fz(ii)*ut2_v;
    ut3=(Fx(ii)*cos(Azi)+Fy(ii)*sin(Azi)).*ut3_h+Fz(ii)*ut3_v;
    
    TT_r(:,:,ii)=ut1;
    TT_t(:,:,ii)=ut2;
    TT_z(:,:,ii)=ut3;
    % convert cylindrical coordinates to Cartesian coordinates
    ux=ut1.*cos(Azi)-ut2.*sin(Azi); %% projection on x-axis
    uy=ut1.*sin(Azi)+ut2.*cos(Azi); %% projection on y-axis
    uz=ut3;
    % stack all sources responses to generate the passive record
    traces_temp=zeros(length(t_total),N_re);
    traces_temp(srt22(ii):srt22(ii)+N-N_shift-N2-1,:)=uz;
    traces_z=traces_z+traces_temp;
    
    traces_temp=zeros(length(t_total),N_re);
    traces_temp(srt22(ii):srt22(ii)+N-N_shift-N2-1,:)=ux;
    traces_x=traces_x+traces_temp;
    
    traces_temp=zeros(length(t_total),N_re);
    traces_temp(srt22(ii):srt22(ii)+N-N_shift-N2-1,:)=uy;
    traces_y=traces_y+traces_temp;
    
hbar.iterate(1); % update progress by one iteration
end
close(hbar); % close the progress bar

traces.x=traces_x;
traces.y=traces_y;
traces.z=traces_z;
traces.t=t_total;
TT.r=TT_r;
TT.t=TT_t;
TT.z=TT_z;
TT.time=t_one;
TT.offset=offset;
TT.fm=fm;
end
